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ABSTRACT 

We present results from axisymmetric, time- dependent hydrodynamical (HD) 
and magnetohydrodynamical (MHD) simulations of a gaseous envelope collapsing 
onto a black hole (BH). We consider gas with so small angular momentum that 
after an initial transient, the flow in the HD case, accretes directly onto a BH 
without forming a rotationally support torus. However, in the MHD case even 
with a very weak initial magnetic field, the flow settles into a configuration with 
four components: (i) an equatorial inflow, (ii) a bipolar outflow, (iii) polar funnel 
outflow, and (iv) polar funnel inflow. We focus our analysis on the second flow 
component of the MHD flow which represents a simple yet robust example of 
a well-organized inflow/outflow solution to the problem of MHD jet formation. 
The jet is heavy, highly magnetized, and driven by magnetic and centrifugal 
forces. A significant fraction of the total energy in the jet is carried out by a 
large scale magnetic field. We review previous simulations, where specific angular 
momentum was higher than that assumed here, and conclude that our bipolar 
outflow develops for a wide range of the properties of the flow near the equator 
and near the poles. Future work on such a simple inflow/outflow solution will 
help to pinpoint the key elements of real jets/outflows as well as help to interpret 
much more complex simulations aimed at studying jet formation and collapse of 
magnetized envelopes. 

Subject headings: accretion, accretion disks - methods: numerical - MHD - stars: 
winds, outflows 



1. Introduction 

The common occurrence and importance of astrophysical jets have stimulated many 
theoretical studies. Both analytic and numerical studies provide strong support for the 
scenario where jets are magnetohydrodynamic (MHD) in character and are associated with 
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an accretion disk, a rotating accretor, or both (e.g., Blandford & Payne 1982; Blandford 1990; 
De Villier et al. 2003). In fact, most MHD simulations of accretion flows show outflows (e.g., 
Uchida & Shibata 1985; Stone & Norman 1994; Hawley & Balbus 2002; De Villier et al. 2004; 
Mizuno et al. 2004; Kato et al. 2004; McKinney & Gammie 2004). 

Here, we describe a study of the time evolution of MHD flows in the vicinity of a 
stationary black hole (BH). This study has been motivated by the results from our previous 
work (Proga & Begelman 2003b, hereafter PB03b; Proga et al 2003) which showed that large- 
scale magnetic fields can produce a jet from a rotationally supported accretion disk or torus 
and also from extremely low angular momentum gas that almost radially accretes onto a BH. 
This falling gas and associated outflow can be well-organized and a simple, self-consistent 
solution for the MHD jet problem. However, in our previous simulations and likely in some 
other simulations, this simple inflow/outflow is a component of a complex convolution of a 
rotationally supported, MHD turbulent torus, the torus corona and outflow. To articulate 
the basic physics that occurs in jet production, we focus on a flow with angular momentum 
so low that, if not for the effects of MHD, the flow would accrete directly onto a BH without 
forming a disk. 



2. Method 

To calculate the flow structure and evolution, we solve the equations of ideal MHD: 

^ + pV.v = 0, (1) 

Dw 1 
p— = -VP - pV$ + — (V x B) x B, (2) 

— =Vx(vxB), (4) 

where p is the mass density, P is the total gas pressure plus radiation pressure, v is the 
fluid velocity, e is the internal energy density, $ is the gravitational potential, and B is the 
magnetic field vector. We adopt an adiabatic equation of state P — (7 — l)e and consider 
models with 7 = 5/3. 
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We perform simulations using the pseudo-Newtonian potential of the central mass 
Qpw — GM/{r — Rs), where R$ = 2GM/c 2 is the Schwarzschild radius, introduced by 
Paczynski & Wiita (1980). This potential approximates general relativistic effects in the 
inner regions, for a nonrotating black hole. In particular, the Paczyhski-Wiita potential 
reproduces the last stable circular orbit at r = 3Rs as well as the marginally bound orbit at 
r = 2R S . 

Our calculations are performed in spherical polar coordinates (r, 9, <fi). We assume axial 
symmetry about the rotational axis of the accretion flow (9 = 0° and 180°). The computa- 
tional domain occupies the radial range = 1.5 Rs < r < r Q = 1200 Rs, and the angular 
range 0° < 9 < 180°. The r — 9 domain is discretized on a non-uniform grid as in PB03b. 

Our calculations use the ZEUS-2D code described by Stone & Norman (1992a,b). We 
adopt PB03b's boundary and initial conditions but we make two modifications to the initial 
conditions (see below). An important element of the initial conditions is that the rotating 
gas has constant specific angular momentum, /, and is confined in a wedge near the equator 
between 9 = 90° — 9 and 9 = 90° + 9 . The wedge of the rotating gas is truncated at 
the radius, r(v r = c^) defined as the radius where the initial radial velocity (as predicted 
by modified Bondi velocity law, see Proga & Begelman 2003a, hereafter PB03a) equals the 
sound speed at infinity, Cqo. 

The two modifications of the PB03b initial conditions are: (i) we set the initial conditions 
exactly as in PB03b and then reduce the density and internal energy by a factor of 100 in the 
part of the computational domain where = (outside the rotating wedge) and (ii) we con- 
sider two field geometries: a purely radial magnetic field as in PB03b and a vertical magnetic 
field defined by the potential A = (A r = 0, Ag = 0, = Ar sin 9). For r sin 9 > r(v r = Cqo), 
we scale the magnitude of the magnetic field using a parameter, f3 Q = ShPb{t ) / B 2 defined 
as the plasma parameter (3 = SttP/B 2 at r OJ so that A = (27rPe(r ) / ' (3 ) 1 ^ 2 (where Pb is the 
gas pressure associated with the Bondi solution at r Q ). For rsin# < r{y r = Cqo), we set the 
constant A to a very small value. 

We introduce these two modifications to reduce the flow complexity. In particular, we 
want the rotating flow to be dominant and we want to follow the evolution of the rotating 
flow with as little interference as possible from other flow components. 

We choose the following units: r = R s , t = 4ir(R 3 s /2GM) 1 / 2 , v = c, B = (Airp^c 2 ) 1 / 2 , 
Po — Poo (the density at infinity for a classic Bondi flow). The force is in units of fo = c 4 /4GM 
and the specific angular momentum is in units of l = 2R s c. 
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3. Results 

To simulate a simple inflow/outflow solution we must appropriately set the model pa- 
rameters. Our initial conditions help to promote the dominance the rotating flow and reduce 
effects of the non-rotating flow. The importance of the non-rotating flow depends not only 
on the initial conditions but also on the minimum value of the density allowed in the simu- 
lations, the so-called density floor, pf. We set pf = (1/r) 1 / 2 . For comparison, the maximum 
density at small radii achieved during the simulation is about 10 4 . To aid the rotating flow in 
reaching an organized solution we set / smaller than the critical angular momentum, 2Rsc. 
[Otherwise the rotating gas will form a rotationally supported torus which becomes turbulent 
and generates a magnetized complex corona and outflow (e.g., PB03b).] Here we describe 
results from the simulations with / = 0.8. We set 9q = 56° (as in a fiducial model in PB03b). 

For our choice of /, the flow near the equator is sub-Keplerian at all radii and will likely 
stay sub-Keplerian during the evolution with appreciable radial velocity (i.e., the equatorial 
flow may continue to accrete supersonically) . This feature distinguishes our simulations from 
many previous ones where the focus was on the outflows from rotationally supported disks 
or torii. Our main focus is on formation of an outflow from nearly radially falling gas away 
from the equator. 

We have performed numerous simulations varying the magnitude and geometry of the 
magnetic field, the specific angular momentum, the numerical resolution and the density 
floor. Here we present results from two models: model A without a magnetic field and 
model B for which (3 = 10 3 and the initial magnetic field is vertical. Model A is a reference 
model which illustrates the flow pattern when accretion proceeds directly whereas model B 
illustrates the dynamics and properties of simplest accretion flow which generates outflow 
without formation of rotationally supported disk and without development of magnetorota- 
tional instability. 

In the early phase of the evolution, when the flow relaxes from the initial conditions, 
both models show the same behaviour. First, the gas near the equator falls in nearly radially 
onto a BH. However, the rotating gas closer to the poles diverges away from the equator 
because of a lack of pressure equilibrium in the 9 direction due to the density and internal 
energy difference between the equatorial wedge and the polar region. The flow at small radii 
becomes gradually radial, with the density decreasing between the equator and the poles. 
After this early phase, the evolution of the two models proceeds differently. In particular, 
the flow in model A settles down into a steady state of direct accretion whereas in model B 
an outflow as well as a direct accretion flow form. 

Figure 1 presents the flow pattern for model A at small radii at t — 3550. The figure 
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shows the density map overplotted with the direction of the poloidal velocity. The flow 
has two components: (i) an equatorial inflow and (ii) a polar funnel inflow. The equatorial 
inflow has non-zero I but its cicularization radius is inside the last stable orbit therefore direct 
accretion occurs. The polar funnel inflow has zero /. The density contrast between the two 
components is due to the fact that at the outer boundary and for the initial conditions the 
density of the non rotating gas is smaller than the density of the rotating gas by a factor 
of 100. At t = 500, the mass accretion rate, M a = 0.08 in units of the Bondi rate (M B ) 
and then gradually increases. At the end of the simulation, t = 3550, M a = 0.123 and still 
continues to grow but very slowly. We estimate that M a will saturate at the level of « 0.125. 
We stress that there is no indication of an outflow in model A. This contrasts with the MHD 
counterpart of this model. 

In model B, the establishment of the equatorial inflow is soon (at t > 715) followed 
by a development of a bipolar outflow from the 'shoulders' of the inflowing gas. At first, 
the outflow is confined to a very narrow range of 9 but with time this range increases. In 
particular, at t=5700, the outflow is at 26° < 6 < 48° and 132° < Q < 154° at r = 20. In 
other words, we observe that with time the boundary between the equatorial inflow and the 
bipolar outflow moves toward the equator. The mass loss rate associated with the outflow, 
Mw increases with time. For example, at t = 1500, Mw = 0.0008 while at t = 5700, 
M w = 0.008 (M w was measured at r = 175). The late time evolution of M a for model B 
differs from that for model A. Namely, for model B, M a = 0.08 at t = 500 and then gradually 
increases to 0.11 at t = 1500. The mass accretion rate stays at this level until t = 5000 and 
then decreases to 0.102 at t = 5700. 

The magnetic field evolves from vertical to nearly radial. This change in the field 
geometry is a natural consequence of radially falling gas because the field lines are dragged 
in with the gas. However, even at the end of the simulation, the field did not evolve into a 
split monopole configuration due to the initial nonradial evolution of the flow. 

Of course, the evolution of the field configuration is accompanied by the growth of the 
field strength. In particular, shear generates the toroidal field in the rotating flow. The 
toroidal field is fastest growing component of the field in the rotating gas and eventually 
dominates the poloidal field, B p = (B r + Bq) 1 ^ 2 . Additionally, the continuously growing 
field becomes dynamically important (we find that an outflow forms in the regions where 
= 8ttP/B^ < 1). In the polar region, the poloidal field is dominant, [3 P = 8nP/B 2 < 1) 
whereas (3$ > 1). 

Figure 2 shows the flow pattern at small radii at t — 5700 for model B. The left and 
right panels show density and \B,p\ maps, respectively. The figure shows also the direction 
of the poloidal velocity and poloidal magnetic field as well as an example of a streamline. 
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The flow has four components: (i) a radial equatorial inflow, (ii) a bipolar outflow (the 
streamline shown in Fig. 2 is typical for this flow component), (iii) a polar funnel outflow, 
and (iv) a polar funnel inflow. The third and fourth components (we call collectively the 
funnel flow) are separated from the inflow/outflow solution by the centrifugal barrier. We 
note that the funnel outflow can be hardly distinguished from the bipolar outflow based on 
the direction of the poloidal velocity. Comparison of the panels of Fig. 2 reveals that across 
the boundary between the bipolar outflow and the funnel flow, three significant changes 
occur: (i) the density dramatically decreases; (ii) B p changes suddenly direction, and (iii) 
increases. We also note that the polar funnel inflow is confined to a very narrow region 
along the rotational axis and its properties/dynamics are difficult to capture by our numerical 
approach because of the imposed density floor and gas heating by artificial viscosity. These 
effects can occasionally lead to production of a thermal outflow in the funnel (e.g., see the 
lower half of the left panel for r < — 11 near the rotational axis). In the remaining part of 
the paper, we focus on the second component of the MHD flow - the inflow/outflow solution. 

The relative simplicity and slow time evolution of the inflow/ outflow solution helps us 
to identify the forces responsible for production of the outflow from the infalling gas. For 
example, we have analyzed and compared all the terms in the equation of motion along 
various streamlines. The important force components in the radial direction are: gravity, f g , 
gas pressure gradient, f p , centrifugal force, f c and gradient of the toroidal field pressure, f m 
(the other components of the Lorentz force are negligible). In the latitudinal direction f p , 
f c and f m determine the total force, f t . 

Figure 3 presents various fluid and magnetic field properties and the forces acting on the 
flow as a function of the pathlength (the pathlength is measured from r = 14 and 9 = 53°) 
along the typical streamline for the inflow/outflow solution. To generate this plot, we assume 
that the solution is steady. This assumption is justified for our purposes here because the 
time changes in the flow occur on the scales significantly longer than the time needed for the 
flow element to follow the segment of the streamline we analyze. 

Fig. 3a shows that the motion in the radial direction is determined in the following 
way: (1) the infall slows down and eventually stops due to the centrifugal force. Up to the 
stagnation point, the fluid specific angular momentum stays nearly constant (dash-dotted 
line in Fig. 3f). The centrifugal force can balance gravity, despite sub-Keplerian rotation on 
the equator, because f c increases from the equator toward the poles [f c oc 1/ (r 3 sin 2 6) for 
/ = const]. Near the stagnation point, two significant changes occur in the force balance: 
(i) f m becomes significant and is directed outward and (ii) f c increases compared to gravity 
because the magnetic torque increases somewhat specific angular momentum of the fluid. 
The combined centrifugal and pressure forces overcome the combined gravity and gas 
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pressure (the latter is directed inwards) and accelerate the flow outward. Although f c > f m 
in the outflow, f c alone is not strong enough to accelerate the flow. Therefore, the flow is 
not 'flung out' by the centrifugal force as in Blandford & Payne's (1982) magnetocentrifugal 
wind but is much more gradually pushed by the pressure of the toroidal field (see Spruit 
1996 and also below for a discussion of heavy loaded MHD winds). 

Fig. 3b shows that as the infalling gas approaches the stagnation point, f p pushes the 
gas away from the equator against the centrifugal force. Near the stagnation point, two 
significant changes occur in the force balance in the 8 direction: (i) f m chances sign and (ii) 
f p increases so that it alone can 'collimate' the gas during the early phase of the outflow 
(this is done against not only f c and also against f m ). 

Farther downstream from the stagnation point, f m again changes sign and becomes 
stronger than f p . Thus, f m appears to play an important role in the outflow collimation at 
large radii whereas f p pushes the outflow away from the equator during the initial accelera- 
tion. However, when the 9 component of f m changes sign for the second time, B r changes 
sign too (the dashed vertical lines in all panels of Fig. 3 mark the location where B r changes 
sign downstream from the stagnation point). In ideal MHD, the magnetic field should be 
dragged in with the flow and no changes in the field orientation should occur within an 
organized outflow due to field freezing. However, in our solution of MHD equations using 
the grid based code, the change of B r sign is due to annihilation (reconnection) of oppositely 
oriented fluxes (i.e., the polar funnel flux vs. the bipolar outflow flux) on the grid-spacing 
scale. We note this numerical limitation and focus on the properties of the bipolar outflow 
upstream from the point where the polar funnel flux annihilates the bipolar outflow flux. 

The bipolar outflow is unbound and it is heavy loaded. To quantify the latter, we 
compute the so-called mass loading parameter /i = {vpV^/v 2 ^) 1 ^ 2 (where va p = {B 2 / Anp) 1 ^ 2 
is the poloidal Alfven spead) introduced by Spruit (1996, see also eq. 2 in Anderson et al. 
2004; note a difference between the definiation of fi in the two references). At the stagnation 
point, /i = 30. Fig. 3d shows that the toroidal field is higher than the poloidal field (compare 
solid, dotted, and dashed lines) while Fig. 3f shows that the outflow does not corotate with 
its base and the outflow does not gain much of specific angular momemtum (see dashed 
and dash-dotted lines in Fig. 3f for the angular velocity, Q and specific angular momentum, 
respectively). All these flow characteristics are consistent with the high value of ji at the 
stagnation point (e.g., Spruit 1996). 

Fig. 3e shows that the flow is super Alfvenic (dotted line). This figure also shows that 
the ratio between the Poynting flux and kinetic energy flux, Fp/Fx increases downstream 
from the stagnation point. The energy flux ratio becomes of order of unity near the location 
of the bipolar outflow 'merging' with the funnel flow (i.e., the location mark by the vertical 
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dashed line). 

4. Summary and Concluding Remarks 

We present results from asymmetric HD and MHD simulations of a rotating gas collaps- 
ing onto a stationary BH. We consider extreme cases where the gas rotation is so small that 
the circularization radius is inside the last stable orbit. Therefore one expects direct accre- 
tion of all gas without formation of a rotationally supported torus no mention of formation 
of an outflow. After an initial transient, the HD flow, as expected, settles into a configura- 
tion where direct accretion occurs and there is not outflow. This contrasts with the MHD 
flow, which even for a very weak initial magnetic field, settles into a configuration with four 
components: (i) a radial equatorial inflow, (ii) a bipolar outflow, (iii) polar funnel outflow, 
and (iv) polar funnel inflow. Our focus here is on the second, unexpected component of the 
MHD flow. 

The bipolar outflow is driven by the magnetic and centrifugal forces. Following the 
streamlines of the bipolar outflow, we identify four stages of the flow motion: (i) a gravi- 
tational collapse, (ii) a gradual slow down and eventual termination of the collapse by the 
centrifugal force, (iii) launching of the outflow by gas pressure which redirects the flow away 
from the equator and by the centrifugal force which redirects the flow away from a BH in the 
radial direction (the latter force is enhanced owning to a magnetic torque spinning up the 
gas), and (iv) an acceleration by the centrifugal force and gradient of the toroidal field pres- 
sure. The outflow is highly magnetized (f3 < 1) and the Poynting flux carries a significant 
fraction of the total energy. 

Despite performing the simulation for thousands of orbits at the inner radius the flow 
did not reach a steady state. One indication of the time evolution is that an increasingly 
larger fraction of the equatorial inflow turns into the bipolar outflow. Additionally, the inflow 
very close to the equator can start turning back, the equatorial symmetry breaks, and the 
flow near the equator starts to circulate. However, not all flow properties evolve with time. 
For example, the angular velocity profile along the equator is time steady and closely follows 
the 1/r scaling as expected for the constant I flow. Generally, the flow evolves over a very 
long time scale into a state where the bipolar outflow 'sandwiches' a smooth radial inflow 
or a complex flow near the equator (our test runs with I = 0.5 show that the equatorial 
flow can remain a smooth inflow for as long as the simulations). We have observed a similar 
bipolar outflow in almost all of our simulations performed recently or performed in PB03b 
and Proga et al. (2003). Our inflow/outflow solution resembles inflow/outflow self-similar 
solutions studied in the context of the core collapse leading to star formation (e.g., see the 
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circulation region in Fig. 1 of Lery el. al 2002; see also Henriksen, R.N., & Valls-Gabaud 
1995, and reference therein). We also note that our solution could be qualitatively similar to 
inflow-outflow circulation found in the simulations of the magnetized cloud contracting under 
self-gravity (Tomisaka 1998) and to the funnel-wall jet found in fully relativistic numerical 
simulations of accretion disks in the Kerr metric (De Villiers t al. 2004). 

We finish with two remarks: (i) The inflow/outflow solution is relevant to several astro- 
physical problems where some of the magnetized collapsing gas may have very low angular 
momentum (e.g., star forming regions, the Galactic center and other supermassive BH envi- 
ronments where the BH is surrounded by the diffused gas and mass losing stars, some GBRs 
and supernovae where it is believed that the explosion is related to a collapse of a massive 
rotating stellar envelope), (ii) Although our physics, initial and boundary conditions are 
simple compared to those in real systems, our inflow/outflow solution is simple and probable 
robust, and will help to pinpoint the key elements of real jets/outflows from low / accretion 
flows. 
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Fig. 1. — A map of logarithmic density overplotted with the direction of the poloidal velocity 
for model A (a zero magnetic field case). 
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Fig. 2. — Maps of logarithmic density (left panel) and toroidal magnetic field (right panel) 
overplotted with an example of a streamline corresponding to an inflow/outflow for model B 
(a non-zero magnetic field case). The maps are also overplotted with the direction of the 
poloidal velocity and the direction of the poloidal field (the left and right panels, respectively). 
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Fig. 3. — From top to bottom: The dominant forces in the radial and latitudinal directions (a and 
b panels, respectively) along the streamline presented in Fig. 2. The fluid and magnetic properties 
as a function of the pathlength along the streamline (c-f panels). For details see the labels and 
the main text. Note that the specific angular momentum, I is reduced by a factor of 10 while the 
radius along the stream line is reduced by a factor of 50 (dot-dashed line and solid line in panel f, 
respectively). The solid vertical lines in all panels correspond to the location where the flow reaches 
the stagnation point (the minimum radius along the streamline). The dashed vertical lines in all 
panels correspond to the location where the radial component of the magnetic field changes sign 
and indicates the location where the magnetic field associated with the bipolar outflow is advected 
into the grid cell with a much stronger, oppositely oriented magnetic field associated with the polar 
funnel. 



